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Abstract 

Several schemes for introducing an artificial dissipation into a central difference approx- 
imation to the Euler and Navier Stokes equations are considered. The focus of the paper is 
on the convective upwind and split pressure (CUSP) scheme, which is designed to support 
single interior point discrete shock waves. This scheme is analyzed and compared in detail 
with scalar dissipation and matrix dissipation ( MATD ) schemes. Resolution capability is 
determined by solving subsonic, transonic, and hypersonic flow problems. A finite-volume 
discretization and a multistage time-stepping scheme with multigrid are used to compute 
solutions to the flow equations. Numerical solutions are also compared with either theoret- 
ical solutions or experimental data. For transonic airfoil flows the best accuracy on coarse 
meshes for aerodynamic coefficients is obtained with a simple MATD scheme. The coarse- 
grid accuracy for the original CUSP scheme is improved by modifying the limiter function 
used with the scheme, giving comparable accuracy to that obtained with the MATD scheme. 
The modifications reduce the background dissipation and provide control over the regions 
where the scheme can become first order. 
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1 Introduction 


Accuracy must be a primary consideration in the construction of any numerical scheme. In 
principle one would like to devise a discrete scheme with the minimum amount of artificial 
dissipation required for stability, as well as convergence in the case of a stationary solution. 
This usually means imposing the additional constraint that the order of the numerical dis- 
sipation is at least one order of magnitude smaller than the desired order of approximation. 
For general fluid dynamic computations the numerical scheme should be designed to have 
high accuracy in smooth regions of the flow held and high resolution at shock waves and con- 
tact discontinuities. According to Harten [3] such discrete formulations, where the accuracy 
away from discontinuities is at least second order, are called high resolution schemes. The 
design of these schemes for systems of conservation laws is generally based on theory devel- 
oped for a scalar conservation law. As a consequence one cannot ensure that the properties 
of the scheme for the scalar equation are valid for the system. In addition, schemes that 
permit high definition of shock waves without oscillations are first order in the neighborhood 
of shocks. Concern naturally arises regarding contamination of the solution, especially in 
the case of viscous flows. For these reasons the properties and resolution capability of this 
class of schemes must be determined through numerical applications for a wide range of flow 
conditions. 

High resolution schemes of particular interest for solving the compressible Euler and 
Navier-Stokes equations are those that allow shock capturing with a single interior point. 
In [6] Jameson presents two schemes with this property that are derived from two different 
forms of flux splitting. One scheme is designated a characteristic split formulation, and 
it employs the flux difference splitting and linearization technique of Roe [19]. With this 
scheme the diffusive flux depends on a flux Jacobian matrix. The other scheme is called the 
convective upwind and split pressure (CUSP) scheme. For this scheme the artificial diffusive 
flux vector associated with a given coordinate direction is expressed in terms of changes in 
the state and flux vectors. A somewhat limited number of inviscid and viscous computations 
have been performed to evaluate these schemes (see [6]-[7] and [27]-[28]). 

We shall investigate and analyze the CUSP scheme, with emphasis on the HCUSP version 
which allows a solution with constant total enthalpy for steady flow. We discuss the shock- 
capturing behavior for various choices of the dissipation coefficients. We introduce a simple 
modification of the limiter function, which is generally used with the scheme, to control 
background dissipation, and thus global accuracy. Global accuracy is also improved by 
introducing parameters into the limiter function to augment control over the regions where 
the CUSP scheme can become first order. The CUSP scheme includes a contribution that is 
scaled according to the local velocity. If the velocity vanishes, as it does for viscous flows, and 
there is a high aspect ratio mesh, the dissipation in the streamwise direction (i.e., direction of 
long side of mesh cell) may not be adequate for convergence. A change in the velocity scaling 
factor based on aspect ratio is presented. The resolution capability of the HCUSP scheme 
is evaluated for subsonic, transonic, and hypersonic flow problems. A detailed comparison 
of the scheme with scalar and matrix dissipation schemes is performed. The scalar scheme 
is based on the dissipation model of Jameson, Schmidt, and Turkel [4], 
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2 Dissipation 


A finite- volume approach is applied to discretize the fluid dynamic equations of motion. The 
computational domain is divided into quadrilateral cells, fixed in time, and for each cell the 
governing equations can be nondimensionalized and written in integral form as follows: 


d_ 

m 


f wdxdy + [ ( fdy — gdx 
£7 J d£l 


vw 

Re 


aa 


(. f v dy - g v dx ), 


( 2 . 1 ) 


where 0 is a generic cell (or cell area) with dfl its boundary. In the scaling factor for the 
viscous terms on the right hand side of (2.1), the quantities 7 , M, and Re are the specific 
heat ratio, Mach number, and Reynolds number, respectively, with M and Re defined in 
terms of nominal conditions. Taking w h k as the cell-averaged solution vector, equation (2.1) 
can be written in semi-discrete form as 


d 

dt 


(^j,k^j,k) T Cwj^fc — 


0 , 


( 2 . 2 ) 


where Qj-y is the area of the cell, and C is a spatial discretization operator defined by 


Id = Cc + Cd + Id AD, 


(2.3) 


with the subscripts (7, D , and AD referring to convection, diffusion, and artificial dissipation. 
In order to simplify the description of the dissipation model, we consider the one-dimensional 
Euler equations of gas dynamics. 


2.1 Scalar Dissipation Model 

The scalar dissipation is based on the model introduced by Jameson, Schmidt, and Turkel 
[4], This model defines a switching function based on a blending of the second and fourth 
differences. The term associated with the operator Cad is expressed as 

Cad w j = —(D 2 — D 4 )wj = dj +1 / 2 — dj_i/ 2 . (2-4) 

Then 

D 2 m, = V [(A j+1/2 Ai/ 2) a ] w i • (2-3) 

D'w, = V [(A )+1/2 T + ’ 1/2 )AVA] (2.6) 

where the index j refers to a cell center, and the operators A and V are forward and backward 
difference operators. The variable scaling factor A is defined as 

^j+ 1/2 = 2 ’ ( 2 -^) 

where A is the largest eigenvalue in absolute value (i.e., spectral radius) of the flux Jacobian 
matrix associated with the Euler equations. For example, in the £ and g directions of 
generalized coordinates (^, 7 ), 

^ = \uy v -vx v \+ c^Jx 2 v + y 2 n , 
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X v = \u,X£ —uyt\ + cyjx 2 ^ + y|. 

The coefficients e ^ and e^ 4 ) use the pressure as a sensor for sharp gradients, and they are 
defined as 


$1/2 = K{2) maX (ffi-l, W ffi + 1, ffi+2), 
Pi - 1 - + ffi+i 


+ 2pj + p J+ i 

£j+i/ 2 = max [°> (« (4) - e f+i/ 2 )] > 


(2.8a) 

(2.8b) 

(2.8c) 


where typical values for the constants k ^ and are in the ranges | to | and -T to gh, 
respectively. We shall refer to (2.4) together with (2.8) as the JST scheme and (2.8) alone as 
the JST switch. The switching function v can be interpreted as a limiter, in the sense that it 
activates the second-difference contribution at extrema and switches off the fourth-difference 
term. Moreover, at shock waves the dissipation is first order, and a first-order upwind scheme 
is produced for a scalar equation. In smooth regions of the flow held the dissipation is third 
order. 

Thus, we have two different dissipation mechanisms at work. The switch determines 
which one is active in any given region. For smooth hows, v is small and the dissipation 
terms consists of a linear fourth difference that damps the high frequencies which the central 
difference scheme does not damp. This is useful for achieving a steady state and is not 
always necessary for time dependent problems [9]. In the neighborhood of large gradients 
in the pressure, v becomes large and switches on the second-difference viscosity while si- 
multaneously reducing the fourth-difference dissipation. This is mainly needed to introduce 
an entropy condition to reduce overshoots near discontinuities and choose the correct shock 
relationships. For subsonic steady state flow this can be turned off by choosing k ^ = 0. 

One possible extension of the scaling factor of (2.7) to multidimensions is isotropic. In 
two dimensions, with (£, r/) denoting arbitrary curvilinear coordinates, the scaling factor 
takes the form 

Xj+1/2 ,k = )i,fc + (^)j+i,fc + (X v )j,k + (^r))j+i,fc]- (2-9) 

Such a scaling is generally satisfactory for inviscid flow problems when typical inviscid flow 
meshes (i.e., cell aspect ratio 0(1)) are used. This factor can cause excessive numerical 
dissipation in cases of meshes with high-aspect-ratio cells. Instead, the scaling factor is 
usually defined as 

Xj+1/2 ,k = )i,fc + (^)j+i,fc]5 (2.10) 


where 


{Xc)j,k — (f>j,k{ r ) (Xt)j,k, (2-11) 

<f>3,k{ r ) = 2 C_1 (1 + 4fc)’ 
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r is the ratio X v /X^ } and the exponent ( is defined by 0 < £ < 1. If £ = 1, the isotropic 
form of (2.9) is recovered. If ( = 0, the scaling in a given direction simply depends on 
the eigenvalue associated with that direction. This scaling is sometimes called individual 
eigenvalue scaling (see [14], [21]). The exponent ( is generally taken to be between | and |. 
Thus, this dissipation scaling factor is between the isotropic and individual eigenvalue scaling 
factors. As demonstrated in [13] and [21], this factor produces a significant improvement in 
accuracy relative to the isotropic factor for high-aspect-ratio meshes, and it permits good 
convergence rates with a multigrid method. 

Using the TVD concept, an alternative for the switch of (2.8) that is TVD for a scalar 
equation is introduced in [22], In one dimension this switch is given by 


\Pj + 1 ~ 2 Pi + Pj- 1 1 
J \Pj+i ~Pj\ + \Pj ~Pj-i\ + e 


( 2 . 12 ) 


and choose 121 = b In practice we usually use a weaker form than (2.12), for example, 


= \Pj + 1 ~ 2 P] + Pj- 1 1 

J (1 — lo)Ttvd + 5 


(2.13) 


where 


Vtvd = \pj + 1 - Pj | + \pj ~ Pj- 1 1, 

V = p J+ 1 + 2 pj +pj-!, 

and 0 < io < 1. The TVD switch of (2.12) is recovered when lu <C 1. Typically lu 1/2. 
In [23] this switch allowed the computation of flows with strong shock waves whereas the 
switch of (2.8) did not. 

2.2 Matrix- Valued Dissipation Model (MATD) 

Sharp resolution of shock waves without oscillations can be achieved by closely imitating an 
upwind scheme in the neighborhood of a shock wave. A key feature of upwind schemes is a 
matrix evaluation of the numerical dissipation. With this evaluation the dissipative terms of 
each discrete equation are scaled by the appropriate eigenvalues of the flux Jacobian matrix 
rather than by the spectral radius, as in the JST scheme. A matrix dissipation model can 
easily be constructed by starting with the JST formulation. 

One can show [22] that the necessary modification of the JST scheme to produce a 
matrix dissipation model is the substitution of |A| for the eigenvalue scaling factor A in 
(2.5) and (2.6). Since the Euler equations are a strongly hyperbolic system, the coefficient 
matrix can be diagonalized. Assume QAQ~ X = A (diagonal matrix). Then |A| is defined as 
\A\ = an d |A| = diag( |Ai|, |A 2 |, |A 3 |), where A; are the forward acoustic, backward 

acoustic, and convective eigenvalues. An efficient way of computing |A| times a vector is 
presented in [22], 

In practice one cannot choose Ai,A 2 ,A 3 as the eigenvalues. Near stagnation points A 3 
approaches zero while near sonic lines Ai or A 2 approaches zero. A zero artificial viscosity 
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would create numerical difficulties. Hence, we limit these values as 

|ffi| = </>max(|Ai|, V n p(A)), |A 2 | = </>max(|A 2 |, V n p(A)), (2-14) 

|A 3 | = </>max(|A 3 |, V e p(A)), (2.15) 

where </> is dehned by (2.11), p(A ) is the spectral radius of A, and the linear eigenvalue A 3 
can be limited differently than the nonlinear eigenvalues. The parameters V n and Vi have 
been determined numerically. Typical values are V n = 0.25 and Vn = 0.025. 

2.3 CUSP Scheme 

In the previous sections we have described the use of an artificial viscosity based on either 
a scalar or matrix coefficient. Inspired by earlier work on flux- vector splitting [34] Liou 
and coworkers designed a scheme called Advection Upstream Splitting Method ( AUSM ) 
[10, 11, 35]. This method was later refined for large-scale 3-D viscous computations in [17]. 
AUSM is based on a splitting of the flux function into convective and pressure contributions. 
In some sense, the pressure terms contribute to the acoustic waves while the velocity terms 
contribute to convective waves. Hence, it is reasonable that these flux terms be treated 
differently. Liou thus considers decompositions of the flux vector that are not based on a 
characteristic decomposition but on Mach number scaled contributions of the left and the 
right states to the interface flux. This decomposition has the disadvantage that it is more 
difficult to develop for other sets of equations compared with a characteristic decomposition. 
A similar type scheme called the Convective Upwind Split Pressure (CUSP) scheme was later 
introduced by Jameson [5] and subsequently modified by Tatsumi, Martinelli, and Jameson 
[7, 8, 27, 28]. The CUSP scheme has some advantages over AUSM. First, one can consider 
the scheme as another type of artificial viscosity, since it is dehned as a sum of the central flux 
average plus a dissipative flux. Hence, it can be readily used with a variety of time-stepping 
schemes (e.g., multistage, LU, implicit, etc.). Second, the CUSP formulation can be used 
in a straightforward manner with multistage schemes which do not evaluate the artificial 
dissipation fluxes at every stage, in order to reduce computational work. Hence, we shall 
only describe the CUSP version of this type of scheme. 

2.3.1 Definition of CUSP Scheme 

Previously, we introduced the scalar and matrix- valued viscosities by considering dj+i/ 2 of 
the form 

dj+ 1/ 2 = ^Qj+ 1 / 2 K +1 - wj). (2.16) 

The factor | is introduced so that we get full upwinding when Q J+ 1 / 2 = /. We note that for 
the scheme to be positive, Q must be sufficiently large. For the matrix viscosity we chose 
Q = \A\ (modified near zero eigenvalues) while for the scalar viscosity we chose Q = a (A) I. 
For the CUSP scheme we instead choose d as a linear combination of w and /. In one 
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dimension we consider two choices for the state vector: 


w = (p pu pE) J 

( p \ ( 0 


f = u 


pu I + 

V pE 


P 
\ up 


and 


w h = (p pu pH y 

l p 


f = U 


pu | + 

V pH 


( ° 

P 

V o 


uw + f p 


UW h + f p . 


The first-order accurate CUSP scheme is defined as 


f (3 

rfj + l/2 = + j - MI,) + -(/, + ! - /,) 


( 2 . 17 ) 


The factor c is included so that v is dimensionless. We thus consider only scalar parameters 
instead of a matrix coefficient, but we have two free parameters, v and [3. The scheme is total 
enthalpy preserving if Wh = (p pu pH) T is chosen as the basis. This choice is denoted HCUSP 
by Jameson [7]. By using the arithmetic average, u = |(u J+ i + Uj), and the definition 


ac = vc + [3u . 


One can rearrange (2.17) to obtain 

1 (3/3 

d 3+i/2 = 2 ac ( w j+i -Wj) + ~(f Pj+1 - f Pj ) + 1 ~ u j)- 

Introducing the Roe matrix Arl, we have /r — /l = Arl(wr — wl). This relation is exact 
if Arl is computed from weighted averages of the left and the right states. That is, 

_ a /PrUr + yfpLU L 

\[pr + ^fpZ 

H _ \/PrH r + ^pEH L 
+ \[Pl 

P = VPRPL 

Then the first-order dissipation is 

dj+i /2 = 7 ^(fiA RL + ucI)(w R - uj l ). (2.19) 

We see from this formula that d is a linear function of A. Recall that |A| is a quadratic 
function of A, by the Cayley-Hamilton Theorem. Hence, it is not possible to bound d by 



( 2 . 18 ) 
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Figure 1: Conditions at shock wave; 

|A|. Since |A| is the minimum dissipation needed for a scheme to be positive [26], the CUSP 
scheme cannot be positive. 

Remark. The parameters (3 and v will be defined later in this section. For these pa- 
rameters the CUSP scheme is not positive at least for M < 1/2. The concepts of TVD and 
positivity were introduced primarily for the treatment of discontinuities. Thus, it is not clear 
theoretically if the loss of positivity for subsonic flow bounded away from the sonic line is 
important. For supersonic flow (M > 1) the CUSP scheme is positive. 

Assume that tin®? subscript L denotes the interior point inside the shock zone, R is the 
state downstream of the shock, and the sta.t® LR is subsonic (as depicted in Figure 1). 
Jameson [7] shows that the downstream point with the state R is in equilibrium if 

VC 

fR-fL + —-(lOR-lOL) = 0. (2.20) 

1 + p 

Substituting the Roe matrix for the difference in / into (2.20) we get 

^Arl + 1 ( Wr — W L ) = 0. 

Hence, wr — wr is an eigenvector of Arr, and —(uc)/( 1 + /3) is the corresponding eigenvalue. 
However, the eigenvalues of Arr are known to be A + , A - and u. If A is an eigenvalue of A, 
then using this formula for uc in (2.19) we have 

c/,+ 1/2 = 7 / [—A/ + (3(A rl - A/)] (ior - w L ). 

In order to have a positive diffusion when u > 0, we require that A be negative (i.e., 
— (nc)/(l + (3) = A - ). Thus, 


uc = — ( 1 + [3 ) A . 


( 2 . 21 ) 
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For u < 0 we obtain similarly 


vc=(l-f3)\ + . (2.22) 

So, we have reduced our two free parameters to one free parameter by demanding a one 
point shock profile. More generally, Jameson shows that one obtains a shock profile with 
one interior point if the following two conditions hold: 

1. When the flow is supersonic through the shock then one obtains a totally upwind flux. 

2. The artificial dissipation Q satisfies a generalized eigenvalue problem 

( Ar L — a RA Q RA ) (w R — w A ) = 0 


at the exit from the shock. 


The second condition is satisfied by both the matrix viscosity and the CUSP scheme; however, 
the scalar viscosity does not satisfy the first condition. We again note that the positive 
condition Q > |A| is satisfied by the scalar and matrix viscosities but not by the CUSP 
viscosity for all Mach numbers. 

What remains to be done is to choose suitable functions for (3 and vc which satisfy the 
above requirements. Jameson’s choice for /3, which is based upon the eigenvalues correspond- 
ing to the acoustic waves, is given by 



+ max (0, “1J_) 

if 

0 < M < 1 


p = < 

max (0, ) 

if 

-1 < M < 0 

(2.23) 


sgn (M) 

if 

\M\ > 1. 



The cutoffs, (3 > 0 for u > 0 and (3 < 0 for u < 0, ensure that the pressure terms are 
discretized centrally for small Mach numbers. Shock capturing with one interior point is 
obtained by taking 


vc = < 


u 


-(1 + (3) A" 
{1-(3)X + 

0 


if (3 = 0 

if (3 > 0 and 0 < M < 1 
if (3 < 0 and — 1 < M < 0 
if \M\ > 1. 


The dissipation coefficients are to be computed from Roe-averaged quantities as in (2.18). 
They provide full upwinding for supersonic flow, (3 = sgn(M), v = 0. The choice of vc = |u| 
for (3 = 0 yields a continuous dissipation coefficient in the subsonic region, and it does 
not smear slip lines with |u| close to zero. This makes the CUSP formulation attractive for 
viscous flow calculations with boundary layers. However, viscous flows are usually discretized 
by using cells with large aspect ratio. It is well known that this situation requires larger 
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dissipation scaling in the direction of the long cell sides than given by |u|. We redefine the 
dissipation coefficients in the individual coordinate directions. For the ^-direction we have 


( max(|rt|, her 
— (1 + (3) A" 


vc£ = r < 


(1 — /3)A+ 

0 


if (3 = 0 
if (3 > 0 and 
0 < M < 1 
if (3 < 0 and 
-1 < M < 0 
if \M\ > 1. 


(2.24) 


where r + and r are functions of the spectral radii in the £ and rj directions (A^ and A^), 
and they are defined as follows: 


r = 


A n 

h 


r + = max(r, 1), r 


min(r, 1). 


The dissipation coefficient in the ^-direction is defined correspondingly. 


2.3.2 Simplified Scheme 

Several modifications of the CUSP scheme have been in use so far. Based upon the Wh = 
(p pu pH) T system the dissipation coefficients presented in [7] and [28] are as follows: 


j \M\ if \M\ > e 

l 2 ( e + ~ T ) if \M\ < e, 


(3 


+ max 

-max(0,^ 

sgn(M) 


if 0 < M < 1 
if -1 < M < 0 
if \M\ > 1. 


(2.25) 


This choice does not allow exact shock capturing because (2.21) and (2.22) are not satisfied. 
Furthermore, Roe averaging has been replaced by arithmetic averaging in [7] and A - , A + by 
u — c, u + c, respectively. This simplification saves a few square roots in the coding of the 
dissipative flux. Equation (2.25) is then 


j \M\ if \M\ > e 

l 2 ( e + ~) ^ \M\ < e, 


(3 


max (0, 2 M — 1) 
< min (0, 2 M + 1) 
sgn(M) 


if 0 < M < 1 
if -1 < M < 0 
if \M\ > 1. 


(2.26) 
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3 Higher Order Scheme 


Having determined vc and /3, we see from (2.17) that the scheme is completely defined in 
terms of w and f p . Formula (2.17), as given, is only first-order accurate, as it depends 
only on d J+ 1 / 2 = Wj + \ — Wj, and so the complete artificial viscosity behaves like a second 
difference. The purpose of this section is to combine a first-order accurate CUSP scheme 
with a high-order dissipation. 

Previously, we considered a combination of a low-order and high-order artificial viscosity 
based on the scalar ( .JST ) switch of (2.8). This switch has the disadvantage that one quan- 
tity, the pressure, controls the shock sensor. Moreover, it forces all variables to be treated 
equal, even though some experience sharp changes through the discontinuity while others 
are continuous across the shock. The requirement to choose a particular flow variable for 
a switch can be eliminated. One can instead limit independently each dependent variable 
in each coordinate direction. Such a limiting allows the construction of a strictly upwind 
scheme for the one-dimensional Euler equations rather than just for a scalar equation. 

In [5] Jameson constructed a family of limiter functions based on the function 


R(u } v) = 1 


u — V 

u \ + M + e 5 


(3.1) 


where q is a positive number and e has the dimensions of u. The parameter e << 1, and in this 
work it is taken to be 10 -1 °. Note that R(u } v) ss 0 whenever u and v have the opposite sign. 
Let w be an element of the solution vector for the governing flow equations. Also, note that 
according to our previous theory [22] R(Aw J+3 / 2} Aw J _ 1 / 2 ) } where Aw J+3 / 2 = w J+ 2 — w J+ 1 , 
would be replaced by iy +1 / 2 , where iy +1 / 2 is the maximum of Vj over the nearest neighbors 
and v is given by (2.12). 

In the results section of this paper we will demonstrate that it can be beneficial to further 
control the regions where the limiter is applied. Hence, we generalize (3.1) to 


R(u } v) = 1 — min(e^, e p , 1) 


u — v 


\u + \ v + e 


with 


e 


V 


0 

5 MUL(M - M lmtt 


if M < Mama 
if M > 




0 


if ^ ^ L'limit 

if V ^ L'limit' 


(3.2) 


The control parameters are the contravariant Mach number M and the pressure switch 
as given in (2.13). Thus, the limiter cannot produce a first-order scheme in regions where 
M < M ii m i t or v < uumit- With the introduction of min(e^, e p , 1) in (3.2) the scheme is not 
very sensitive to the value of the exponent q (typically, q = 1 or q = 2). 
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Define the limiter function L(u } v) by 


L(u,v) = R(u,v ) — ^ — . (3-3) 

At the mesh cell interface j + 1/2, we define the left and right states for each dependent 
variable as 

w L = w 3 + l-L(Aw J+3/2 , Aw 3 _ 1/2 ), 

1 (3-4) 

w R = w 3+1 - -L(Aw 3+3/2} Aw 3 _ 1/2 ) } 


and so 


w R - w L = Aw 3+l / 2 - L(Au) i+3 / 2 , A^_ 1/2 ). (3.5) 

For the artificial viscosity all differences will be based on wr — wr. In the neighborhood of 
shock waves R(u } v ) and hence L(u,v ) are close to zero. Moreover, wr — wr = Auy +1 / 2 , 
resulting in a first-order scheme for the artificial viscosity. For smooth flow R(u } v) = 1, and 
L(u } v) = (u + v)/2. Hence, in a smooth region 

WR ~ Wr =Aw j + 1 / 2 ~ L(Aw j+ 3 / 2 ,Aw j _ 1 / 2 ) 

a Aw ]+3/2 + A W 3 _i/2 

~ AW J + 1/2 — ± ~ 

— ^ A 3 

- ~2 Wj+1 / 2 ' 

Thus, in the smooth regions wr — wr behaves as a third difference, while in the vicinity of 
shock waves it behaves as a first difference. Consequently, (3.5) has similar properties to the 
■JST scheme. One can obtain the relationship between (3.5) and the .JST scheme by defining 
the diffusive flux d 3+ i/ 2 as 

dj + 1/2 = a j + l/2 (wr — Wr) , a J + 1/2 = A J + 1 / 2 , (3-7) 

where k ^ is a parameter, and A is the spectral radius of the associated flux Jacobian matrix. 

One difference between the JST scheme and (3.5) involves the parameters k ^ and k ^ 
for the second and fourth differences, respectively. Both k ^ and are free parameters 
in the JST scheme. As seen from (3.6) and (3.7) these parameters are automatically chosen 
as k/ 2 ) and with (3.5). Furthermore, for the matrix viscosity (see Section 2.2) and 

the CUSP scheme (described in Section 2.3) k ^ and so we no longer have any free 

parameters. The coefficient of the second difference is chosen as | so that the scheme is fully 
upwind for supersonic flows. However, the fourth-difference viscosity is introduced only to 
accelerate the convergence to a steady state by eliminating the decoupling of the odd and 
even points. Hence, we wish k ^ to be as small as possible for accuracy while still achieving 
a good convergence rate. It does not seem reasonable to connect the two components of the 
artificial viscosity. In section 4 we will compare the magnitude of the scalar viscosity and 
the CUSP scheme. 


( 3 . 6 ) 
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One can generalize the limiter function of (3.3) by reintroducing a free parameter that 
essentially governs the level of the third-order viscosity in the smooth regions. The resulting 
scheme has the disadvantage that a free parameter must be chosen; however, it has the 
advantages of greater flexibility and increased accuracy. We now define a new L as 

L(u } v } w) = R(u } w) ■ (1 — 4:H^)v + 4/J 4 ) — — , (3-8) 

where left and right state values are determined by 

w L = Wj + ^L(Aw j+3 / 2 , Auy+i/ 2 , Arc J _ 1/2 ), 
w R = w J+ 1 - ^L(Aw j+3 / 2 , Aw j+1/2 , Aw 3 _ 1/2 ). 

When = |, the L of (3.8) reduces to the original L of (3.3). At shock waves R(u } w) ss 0, 
and we again have w R — w R = Aw j+ i/ 2 . For smooth regions of the flow held we have 
wr-w l = -2k^A 3 w j+1 / 2 . 

One difficulty with (3.1), and indeed with any TVD switch, is that it limits the differences 
near minima and maxima independent of the amplitude of the function. Hence, in the far 
held where the solution is almost uniform the low-order scheme is activated by small noise 
levels. Since this occurs randomly it frequently prevents the convergence of the residual 
beyond three or four orders of magnitude. The use of (3.2) eliminates this difficulty. 

The extrapolation technique of (3.8) can be used with (2.17) to get the hrst difference to 
higher order accuracy. Then, the states corresponding to higher order accuracy are obtained 
in a way similar to van Leer’s MUSCL approach [34], To impose monotonicity one can apply 
the limiter discussed in this section. For example, we can replace (2.17) by 

vc j3 

d 3 + 1/2 = Y^ Wr ~ Wl ^ + 2 ~ fp( WL )} ’ 

where w R} wl are given by (3.4). This procedure was followed throughout the numerical 
examples presented in Section 7. Application of (3.4) to the Wh = (/? pu pH) T variables 
still allows total enthalpy to be preserved in the higher order scheme. When a multigrid 
algorithm is used to solve the governing flow equations, the higher order scheme is applied 
only on the finest mesh, and the lower order scheme is applied on the coarser meshes. 
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4 Analysis of CUSP Scheme 

The eigenvalues of (3Arl + vcl (see 2.19) are /.q c, /i 2 c , and /i 3 c. Using the simplifications 
of (2.26) the eigenvalues are: 

hi = 1^1, 

' \M\ if \M\ < \ 

/i 2 = < a + (3 if | < M < 1 

|M + 1| if |M| > 1, 

' \M\ if \M\ < \ 

/i 3 = < a, — (3 if \ A M < 1 

|M — 1 1 if |M| > 1 

We note that for \M\ < | all three eigenvalues of the artihcial viscosity are equal, and so we 
have the equivalent of a scalar viscosity. The scalar viscosity now scales with \M\c rather than 
(\M\ + 1 )c as in the JST scalar viscosity. This is more similar to the case of preconditioning 
where all the eigenvalues are approximately \M\c for low speed flow. Hence, we expect that 
the CUSP dissipation should work properly for very low Mach numbers provided the central 
flux terms are augmented by a suitable preconditioning matrix. 

In the subsonic range where (3 = 0, all of the versions of the CUSP scheme do not 
satisfy (2.21) and (2.22), which are necessary for shock capturing. Thus the cell-face Mach 
numbers in the shock structure have to be larger than about 0.5 in order to avoid post-shock 
oscillations. The motivation to design uc = |rt| when (3 = 0 has already been discussed. 
However, the choice of the function for (3 , as given in (2.23), is not necessarily optimal. For 
example, (3 = max( 0, ( u + ~)/(u — A - )) would allow shock capturing for Mach numbers 

down to about 1/3, but the subsonic dissipation would be twice as large, uc = 2|rt| for (3 = 0. 
Nevertheless, our own experience gained from a number of numerical applications suggests 
that there is no need for further modifications of (3. 

It is rather difficult to compare the effect of the parameter k/ 4 ) of the JST and the CUSP 
schemes, since these schemes also include eigenvalue information which is not the same in the 
two schemes. To isolate the effect we consider a low Mach number flow with preconditioning 
(see [31] for details). Now both switches are based on the convective eigenvalue u. A typical 
value for the JST scheme is k( 4 ) = A However, for an aspect ratio of one the Martinelli 
scaling [13] adds another factor of two. The parameter a = 0 in the preconditioning adds an 
additional factor of approximately 2.6. Hence, the effective constant multiplying the fourth 
difference is about which is somewhat smaller than the | used with the original CUSP 
scheme. For transonic flows it is more difficult to compare the levels of dissipation. However, 
it seems that the original CUSP scheme yields too high a viscosity level and so the k ( 4 ) 
introduced in (3.8) should be reduced to less than Numerical computations demonstrate 
the improved accuracy (though slower convergence) for standard transonic turbulent flows 
when k W is reduced. 


14 



5 Low Speed Preconditioning 

For low Mach numbers standard algorithms converge to a steady state very slowly because of 
the disparity between the convective and acoustic eigenvalues. Furthermore, it is also found 
that most schemes give very poor results for these low Mach number flows, even when a 
steady state is achieved ([30]- [32]). One way to overcome these difficulties is to precondition 
the equations by multiplying the time derivative terms by a matrix P _1 . If P is appropriately 
chosen, then one can reduce the disparity of the wave speeds. In this section preconditioning 
is applied to the different dissipation schemes. Details of preconditioning techniques are 
given in [31, 30]. 

To simplify the presentation we shall only consider a one-dimensional system. The ex- 
tension to multidimensions is straightforward. Consider the system of equations 

dw df 

b — = 0. 

c)t dx 

We replace this by the preconditioned system 


P~ 1 — + ^ L = 0 


dw df 


or in quasi-linear form 


dt 

,-i dw 


dx 


, dw 

m + A ir x = 0 


where A = df /dw. Introducing an artificial viscosity in conservation form, we get 

+ Af = 2-V b l2 »P-U(P.4)A„] = < t+ 1 /2 - 4-1/2 

dt 2Ax Ax L J Ax 


(5.1) 


dw p _Do / 
dt 2Ax 


= fv [e (2) P _1 F(PA)An;] = P ^' +1/2 dj 1/2 


where D 0 denotes a central difference, V is a backward difference, and A is a forward 
difference operator. 

We first consider the matrix- valued viscosity, and thus F(PA) = |PA|. The artificial 
viscosity is 

dj+ 1/2 = e y+i /2 -^j-i-i/ 2 1 (-E >j/ 4-)i+i/2 1 (^i+i — w j) 

Pjdj+1/2 = — w i) m 

When PA has only three distinct eigenvalues, then by the Cayley-Hamilton theorem |PA| = 
a 0 I + ceiPA A a 2 (PA) 2 } where the coefficients cq depend on the eigenvalues of PA. So 

P -1 |PA| = a 0 P _1 + cciA + a 2 PA 2 


We next consider the CUSP artificial viscosity. The artificial viscosity term is given by 
a 0 Aw A ol\ Af ~ a 0 Aw A aiAAw with the appropriate coefficients cq for the CUSP scheme. 


15 



This has the same form as the matrix-valued artificial viscosity without the quadratic term, 
and so by the identical reasoning we get 


1 fj 

d j +1/2 = -ucP j l 1 / 2 (w j+1 - Wj ) + -(fj+i - fj) 


P ,d 


(5.2) 


P, 


3^3 


+1/2 - ^ i/c PjP j +i/ 2 (^j+i _ w j) + trPj(/j+i _ fj) 


(compare with (2.17)). In theory the parameters v and (3 should depend on the eigenvalues 
of PA rather than A and so are not the same as the nonpreconditioned version of CUSP. 
However, in order to not interfere with the shock properties of the CUSP scheme we turn 
off the preconditioning for M > 1/2. Hence, the relationship (2.21) is still valid. The 
parameter (3 would then be chosen by (2.23) where the eigenvalues A + ,A“ should account 
for the preconditioning. In addition, for \M\ < 1/2 the CUSP scheme reduces to a scalar 
viscosity proportional to the convective velocity which is appropriate for preconditioning. 
Hence, it is reasonable to use the same parameters v and (3 for the preconditioned CUSP 
as given by (2.21) based on the original eigenvalues or one of the simplifications previously 
discussed. The advantages of combining the CUSP scheme with preconditioning are shown 
in Section 7. Additional results with the preconditioned CUSP scheme are presented in [33]. 


6 Summary 

The central difference scheme requires an artificial viscosity in order to both prevent oscil- 
lations near shocks and damp high frequencies, enabling the iteration procedure to reach a 
steady state. In the Jameson, Schmidt, Turkel ( .JST ) formulation these artificial viscosities 
are provided by second and fourth differences of the variables with a scalar coefficient in- 
cluded. This scalar coefficient depends on the largest eigenvalue (in each direction) to scale 
the size of the viscosity. In addition, the coefficient depends on the second difference of the 
pressure to sense shocks. In the neighborhood of shocks the fourth difference is turned off 
while the second difference prevents overshoots. In smooth regions of the flow the second 
difference (which leads to first-order accuracy) is minimal while the fourth difference damps 
the high-frequency errors. 

This technique works quite well for transonic flow and was the main approach for many 
years. With the increasing popularity of upwind schemes it was seen that this scheme is 
less accurate than upwind schemes, especially on coarse meshes (see e.g., [1]). This led 
to the introduction of a matrix-valued coefficient in the artificial viscosity (dissipation) that 
mimics the effects of an upwind scheme, but within the context of a central difference scheme 
with an artificial dissipation coupled with a multistage time advancement. Later Jameson 
introduced the CUSP scheme, which is in between the matrix dissipation and the scalar 
dissipation schemes. With the CUSP scheme the dissipation is a function of the Mach 
number and becomes fully upwind in supersonic regions similar to the matrix dissipation. 
However, it avoids the need for a full-matrix coefficient while still obtaining one point shock 
profiles. The CUSP scheme is more expensive than the matrix-viscosity method, since it 
uses extrapolation for each flow variable and limiting depending on the variable to achieve 
second-order accuracy, 
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The idea of using the Mach number to mimic fully upwind methods has been useful in 
other applications besides the artificial viscosity. One can use the Mach number to adjust the 
parameters of the residual smoothing so that it becomes fully upwind in supersonic regions. 
Similarly one can construct a multigrid method with weighting factors depending on the 
Mach number such that it becomes fully upwind in supersonic regions. We label this 

Poor Man’s Upwinding 


Advantages: 

• Cheap 

• Fully upwind in supersonic flow 

Disadvantage: 

• Not fully upwind for subsonic flow 
Applications 

• Residual smoothing 

• CUSP scheme 

• Multigrid 

In the opposite direction the Mach number can be used to construct a preconditioning 
that is useful in low Mach number regions. Then the Mach number can be used to turn off 
the preconditioning in supersonic regions. 

7 Numerical Results 

In the numerical applications presented here we assess the accuracy and shock capturing 
capabilities of the CUSP scheme. Since the version of the CUSP scheme that we consider is 
expressed in terms of the total enthalpy H and is H preserving for inviscid flows, it is usually 
called the HCUSP scheme. Comparisons are made between the HCUSP and MATD schemes. 
The commonly used scalar dissipation scheme is also included in some of the comparisons. 
In so doing one can clearly see the superiority of the high-resolution HCUSP and MATD 
schemes on even relatively coarse meshes (i.e., 8 cells in the boundary layer of a viscous 
flow). The flow problems considered in the evaluation of these numerical diffusion schemes 
include the following: 1) Inviscid flow over airfoils, 2) laminar flow over a flat plate, 3) 
turbulent flow over an airfoil, 4) inviscid and viscous hypersonic flow over a 2-D wedge. The 
computational effort and convergence behavior in computing these solutions are given. In 
all cases a five-stage Runge-Kutta scheme in conjunction with the convergence acceleration 
techniques of local time stepping, implicit residual smoothing, and multigrid was used. 
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Figure 2a: Inviscid pressure distributions computed with MATD scheme (NAGA 0012 airfoil, 
M = 0.80, a = 1.25°). 



Figure 2b: Inviscid pressure distributions computed with HCUSP scheme (NAGA 0012 
airfoil, M = 0.80, a = 1.25°). 
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x/c 

Figure 3a: MATD solutions near lower surface shock (NAGA 0012 airfoil, M = 0.80, a = 
1.25°). 



X/C 


Figure 3b: HCUSP solutions near lower surface shock (NAGA 0012 airfoil, M = 0.80, 
a = 1.25°). 
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Dissipation 

Scheme 

Grid 

C L 

Cd 

MATD 

192 x 32 

0.3521 

0.02249 


384 x 64 

0.3550 

0.02256 


768 x 128 

0.3552 

0.02256 

HCUSP 

192 x 32 

0.3667 

0.02419 


384 x 64 

0.3610 

0.02310 


768 x 128 

0.3582 

0.02278 

HCUSP 

192 x 32 

0.3639 

0.02297 

(modihed) 

384 x 64 

0.3592 

0.02279 


768 x 128 

0.3563 

0.02269 


Table I: Lift and drag coefficients for inviscid flow over NACA 0012 airfoil ( M ^ = 0.80, 
a = 1.25°). 

The first case is similar to the application published in [7]. Results obtained with the 
MATD and HCUSP schemes for inviscid transonic flow over the NACA 0012 airfoil are 
compared in Figures 2 and 3. The free-stream Mach number for this case is 0.8 and the 
angle of attack is 1.25°. Solutions were computed on three successively finer C-topology 
meshes. The coarsest mesh contained 192 X 32 cells, with 160 cells on the airfoil, and for each 
sequential mesh the number of cells in each coordinate direction was doubled. The principal 
differences between the solutions occur at the shock waves. Since the MATD scheme uses 
a pressure switch for all the flow equations, it cannot capture a shock with a single interior 
point. It requires three interior points. Nevertheless, the resolution of the stronger upper 
surface shock is nearly the same for both the MATD and HCUSP schemes on the 384 X 64 and 
768 X 128 meshes. With the MATD formulation there is some smearing on the 192 X 32 mesh. 
As is clearly evident in Figure 3 the HCUSP scheme allows a sharp definition of the weak 
lower surface shock and the Zierep singularity that immediately follows. The aerodynamic 
coefficients calculated with the MATD , original HCUSP (with limiter of (3.1)), and modified 
HCUSP (with limiter of (3.2)) schemes are presented in Table I. The coefficients computed 
with the MATD scheme on the 384 X 64 mesh essentially agree with those for the finest 
grid. The lift coefficients determined with the original and modified HCUSP schemes on 
the corresponding meshes are slightly higher, with the finest grid values approaching those 
obtained with the MATD scheme. Drag coefficients obtained with the modified HCUSP 
scheme are in closer agreement with those obtained with the MATD scheme, especially on 
the coarsest grid. Later, in the discussion viscous airfoil flow results we will show the behavior 
of the two forms of the limiter in the flow held. 

As an initial evaluation of the dissipation schemes for viscous hows we consider low-speed 
(Mqo = 0.15) how over a hat plate at zero incidence. For this how the Reynolds number per 
unit length is 10 5 . The computational domain is a rectangle. With respect to the leading 
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edge of the plate, the domain extends two plate lengths upstream and one plate length 
downstream. The upper boundary is four plate lengths above the plate. Solutions were 
computed on the same domain and grids used in [28]. Starting with the finest mesh, coarser 
meshes were determined by successively eliminating every other mesh line. The finest grid 
consists of 512 X 128 cells, with 384 cells on the plate. In the direction y normal to the plate 
the grid is spaced uniformly in the boundary-layer coordinate r) (y = y / Re l J 2 ), where x is the 
coordinate parallel to the surface, and Re x is the Reynolds number based on distance from 
the leading edge of the plate). Thus, there is constant resolution of the boundary layer at 
each location along the plate. Outside the boundary layer the grid is stretched exponentially. 
In order to resolve the region in the vicinity of the stagnation point, the grid is clustered 
at the leading edge of the plate. At the surface of the plate no-slip and adiabatic boundary 
conditions are enforced. Along the boundary upstream of the leading edge, a symmetry 
condition is applied. Characteristic type boundary conditions are used at the upstream, 
downstream, and upper boundaries. 

A comparison of the velocity profile at X/ L = 0.82 computed with the scalar, matrix, 
and HCUSP dissipation forms is displayed in Figure 4. Even with just eight points in the 
boundary layer (64 X 16 grid) the MATD and HCUSP schemes nearly replicate the Blasius 
solution. As demonstrated in [1] scalar dissipation can produce serious contamination. With 
the scalar dissipation, more than 32 points are required in the boundary layer to obtain a 
grid converged solution. For the MATD and HCUSP schemes the variation of the errors 
(relative to the Blasius solution) in the calculated skin friction, displacement thickness, and 
momentum thickness are shown in Figures 5a and 5b. The standard definitions given in [20] 
are used for these boundary-layer quantities. The errors in all the boundary-layer parameters 
are quite similar for the high-resolution schemes. This is not surprising since both schemes 
have a scaling factor that vanishes as the surface is approached. 

Transonic flow over the RAE 2822 airfoil is the next test case. The free-stream Mach 
number is 0.73, the angle of attack is 2.79°, and the Reynolds number based on the airfoil 
chord is 6.5 X 10 6 . Transition of the flow from laminar to turbulent is fixed at the 3% chord 
location. The C-type grids used in the computations are as follows: (1) 160 X 32 with 128 
cells on the airfoil, (2) 320 x 64 with 256 cells on the airfoil, and (3) 640 x 128 with 512 
cells on the airfoil. In order to determine the effect of further mesh refinement a calculation 
was performed with the MATD scheme on a 1280 X 256 grid. As in the flat-plate case, 
each successively coarser grid was generated by eliminating every other mesh line in both 
coordinate directions of the finer mesh. The outer boundary is located 20 chords from the 
airfoil. The normal spacing at the surface of the 640 X 128 mesh is 7.5 X 10 -6 chords. At the 
leading and trailing edges of the airfoil the mesh is clustered, giving tangential spacings of 
1.17 X 10 -3 and 1.86 X 10 -3 chords, respectively. These critical mesh-defining spacings are 
roughly doubled with each mesh coarsening. 

In Figure 6 the pressure (C p ) and surface skin-friction (C /) distributions computed with 
the different dissipation schemes for the 160 X 32 mesh described are shown, along with the 
experimental data of [2], As in the inviscid cases the primary differences in the solutions occur 
at the shock wave. Both the scalar dissipation ( SCALAR ) and HCUSP schemes produce 
a solution with the shock too far upstream. This is an unexpected result for the HCUSP 
scheme. The acceleration of the flow upstream of the shock is underpredicted relative to the 
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= Y*SQRT(Re x )/X r\ = Y*SQRT(Re x )/X r| = Y*SQRT(Re x )/X 




4a: Tangential and transverse velocity profiles, X/L = 0.82, 64 X 16 grid. 




4b: Tangential and transverse velocity profiles, X/L = 0.82, 128 X 32 grid. 




4c: Tangential and transverse velocity profiles, X/L = 0.82, 256 X 64 grid. 


Figure 4: Boundary-layer profiles on flat plate with M = 0.15 and Re = 10 5 . 
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Figure 5a: Comparison of results with the MATD scheme and Blasius solution (M = 0.15 
and Re = 10 5 ). 



Figure 5b: Comparison of results with the HCUSP scheme and Blasius solution (M = 0.15 
and Re = 10 5 ). 
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Figure 6a: Comparison of pressure distributions with SCALAR , MATD , and HCUSP 

schemes on 160 x 32 grid (RAE 2822 airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 



Figure 6b: Comparison of skin-friction distributions with SCALAR , MATD , and HCUSP 
schemes on 160 X 32 grid (RAE 2822 airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 
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RAE 2822 airfoil 
grid 160x32 



(b) 


Figure 7: Contours of function R(u, v) in limiter used with HCUSP scheme; (a) basic limiter 
and (b) modified limiter that depends on the contravariant Mach number and pressure 
switch. 
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finest grid. In [24] the adverse effect of a smooth limiter on the accuracy of the solution in 
the vicinity of flow transition, and thus on the acceleration of the flow upstream of the shock, 
is demonstrated. Therefore, such a result with the HCUSP scheme could be a consequence 
of the smooth limiter being used. Thus, we examined the behavior of the limiter in the flow 
held. 

The action of the limiter is revealed by the contour plot of Figure 7 for the minimum of the 
limiter function R(u } v ) (see (3.2)) taken over all four how variables. The contours indicate 
that the basic limiter produces a hrst-order scheme over signihcant portions of the how 
held. This result suggests that the inaccuracy on the coarse grid with the HCUSP scheme 
is not simply a consequence of the behavior of the limiter in the transition region. Figure 
7 also shows contours of the modihed function R(u } v) which uses both the contravariant 
Mach number and the pressure switch of (2.13). With this function the low-order scheme 
occurs only at shock waves. Coarse grid results obtained with the basic and modihed limiter 
functions are displayed in Figure 8. The shock locations computed with the modihed HCUSP 
scheme and the MATD scheme are nearly the same. 

In Figures 9 and 10 the solutions computed on the hner grids with the modihed HCUSP 
scheme are compared with the other dissipation schemes. The pressure and skin-friction 
distributions obtained with the MATD and modihed HCUSP schemes exhibit little difference 
on each mesh. The SCALAR scheme begins to show fairly close agreement with those from 
the other schemes only on the 640 X 128 grid. With both the SCALAR and the MATD 
schemes a nonphysical increase in the skin-friction solution on the upper surface appears 
at the trailing edge of the airfoil. This nonphysical increase is caused primarily by the 
aspect-ratio function of (2.11). As evident in Figure 11, this behavior does not occur in 
the solution obtained with the HCUSP scheme. The computed aerodynamic coefficients, 
including the pressure and friction contributions to the total drag, are given in Table II. 
On each mesh the lift and drag coefficients corresponding to the solution obtained with the 
MATD scheme exhibit the closest agreement with the 1280 X 256 grid values. There are only 
small discrepancies in the coefficients associated with the MATD and the modihed HCUSP 
schemes on the 320 X 64 grid (see also Figure 12). 

Convergence behavior for the HCUSP and MATD schemes is similar. For each scheme five 
levels of multigrid were used and either 50 or 100 cycles were executed on two coarser meshes 
in order to obtain an initial solution. On the 320 X 64 grid the average rate of reduction 
of the residual with both schemes is about 0.92 for 100 cycles on the finest mesh. Figure 
13 shows the effect of modifying the limiter according to (3.8) and (3.2) on the convergence 
with the HCUSP scheme. It also indicates the effect of the modification given by (2.24) to 
vc in the HCUSP scheme. The convergence is improved by using the 2-D formulation for the 
dissipation coefficient vc. Convergence stall can occur with the original limiter. With the 
modihed limiter and the pressure switch this stall is prevented. Note that convergence with 
( = 0 was possible for this transonic case but not for the hypersonic case presented below. 

The fourth case is the hypersonic 2-D how over a blunt wedge. Figure 14 displays the 
second-order accurate solutions obtained for viscous and inviscid how by using identical 
meshes of 64 X 48 cells. Physical diffusion is so large that the shock prohle is signihcantly 
smeared in the viscous result. For inviscid how, on the other hand, we obtain perfect 
capturing with a single interior point in the shock structure by using the formulation of 
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Figure 8a: Effect of modifications (pressure switch and reduced background dissipation) 
in HCUSP scheme on pressure (160 x 32 grid, RAE 2822 airfoil, M = 0.73, a = 2.79°, 
Re = 6.5 x 10 6 ). 



X/C 

Figure 8b: Effect of modifications (pressure switch and reduced background dissipation) in 
HCUSP scheme on skin friction (160 x 32 grid, RAE 2822 airfoil, M = 0.73, a = 2.79°, 
Re = 6.5 x 10 6 ). 
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Figure 9a: Comparison of pressure distributions with SCALAR , MATD , and HCUSP 

schemes on 320 x 64 grid (RAE 2822 airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 



Figure 9b: Comparison of skin-friction distributions with SCALAR , MATD , and HCUSP 
schemes on 320 x 64 grid (RAE 2822 airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 
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Figure 10a: Comparison of pressure distributions with SCALAR , MATD , and HCUSP 
schemes on 640 x 128 grid (RAE 2822 airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 



Figure 10b: Comparison of skin-friction distributions with SCALAR , MATD , and HCUSP 
schemes on 640 x 128 grid (RAE 2822 airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 
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Figure 11: Behavior of skin-friction at airfoil trailing edge with SCALAR , MATD , and 
HCUSP schemes on 320 x 64 grid (RAE 2822 airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 


Dissipation 

Scheme 

Grid 

C L 

C'D 

C Dp 

C'D f 

SCALAR 

160 x 32 

0.8172 

0.01728 

0.01275 

0.004532 


320 x 64 

0.8331 

0.01743 

0.01194 

0.005487 


640 x 128 

0.8532 

0.01782 

0.01225 

0.005574 

MATD 

160 x 32 

0.8304 

0.01818 

0.01251 

0.005662 


320 x 64 

0.8538 

0.01808 

0.01250 

0.005571 


640 x 128 

0.8597 

0.01799 

0.01246 

0.005535 


1280 x 256 

0.8611 

0.01800 

0.01246 

0.005544 

HCUSP 

160 x 32 

0.7987 

0.01926 

0.01367 

0.005594 


320 x 64 

0.8493 

0.01831 

0.01263 

0.005679 


640 x 128 

0.8592 

0.01803 

0.01245 

0.005585 

HCUSP 

160 x 32 

0.8271 

0.01760 

0.01190 

0.005701 

(modified) 

320 x 64 

0.8565 

0.01801 

0.01234 

0.005673 


640 x 128 

0.8604 

0.01798 

0.01240 

0.005581 


Table II: Lift and drag coefficients for turbulent flow over RAE 2822 airfoil {M lZC , = 0.73, 
a = 2.79°, Re c = 6.5 x 10 6 ). 
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Figure 12a: Variation of lift coefficient with reciprocal of number of points (RAE 2822 airfoil, 
M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 



Figure 12b: Variation of drag coefficient with reciprocal of number of points (RAE 2822 
airfoil, M = 0.73, a = 2.79°, Re = 6.5 x 10 6 ). 
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Figure 13a: Original limiter. 



Figure 13b: Modified limiter. 


Figure 13. Effect of limiter and modified uc on convergence history of HCUSP scheme (RAE 
2822 airfoil, \l . = 0.73, a = 2.79°, Re c = 6.5 x 10 6 ). 
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\ M =10, a=0° 

\ \yr>5, Re D =10000 
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(2.23) and (2.24). Detailed comparisons of the hypersonic wedge flow solutions yielded by 
the CUSP scheme and AUSM have been presented in [16]. It was found that the shock 
capturing capabilities of both schemes are essentially equal. A comparison of shock profiles 
for the exact and the simplified coefficients is given in Figure 15. Here, we have chosen 
the first-order scheme in order to address the pure shock capturing capability of the CUSP 
scheme without interference from the limiter. The simplified dissipation coefficients of (2.26) 
produce strong oscillations at the shock, even though there is substantial physical diffusion 
present. Hence, it is concluded that an accurate implementation of dissipation coefficients 
is a requirement for hypersonic flows with strong shocks. 

Some applications of the MATD scheme to hypersonic flow problems are given in [23]. 
However, we fold that matrix dissipation combined with a pressure-based sensor in order to 
switch from second to fourth differences has not yet resulted in sufficient robustness to deal 
with hypersonic flow phenomena in general. In particular, it seems that the user defined 
coefficients in (2.13) - (2.15) need adjustment depending on the flow problem. Moreover, it is 
well known that matrix dissipation schemes suffer from an instability known as the carbuncle 
problem [15], and they need rather large values of V n and V] in order to restore stability. 

The final set of results show the behavior of the HCUSP scheme with preconditioning. 
Inviscid solutions for flow over a NACA 0012 airfoil were computed on a C-type grid with 
224 X 40 cells and clustering at the leading and trailing edges. In Figure 16 Mach number 
contours delineate the effect of the free-stream Mach number on the solutions obtained with 
the preconditioned HCUSP scheme. Figure 17 clearly illustrates the benefits of precondi- 
tioning on the HCUSP scheme. There is a substantial improvement in not only the quality 
of the solution but also the convergence behavior with the scheme. 

Comparisons of computation times indicate that the HCUSP scheme needs about 25% 
more computer time than the basic scalar dissipation of Section 2.1. The MATD scheme only 
requires about 15% additional time. This reduction is primarily a consequence of the single 
evaluation of the limiter function. Due to lower inherent dissipation, computations with the 
HCUSP formulation converge somewhat slower for transonic flows than those with simple 
scalar dissipation. The major advantage of the HCUSP approach is that it is more accurate 
and more robust than scalar viscosity. Our numerical tests indicate that the accuracy of 
the CUSP scheme is close that of matrix dissipation for transonic flows provided the first- 
order scheme is activated at shock waves only. For hypersonic flows it seems to be more 
robust than the matrix viscosity even though it is not positive. Since the HCUSP scheme 
is implemented through artificial dissipative terms, it does not have to be applied at each 
stage of the Runge-Kutta method. In particular, the diffusive fluxes can be evaluated only 
at the first, third, and fifth stages of a five-stage method, as is typically done for the scalar 
dissipation. 
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center line flow 



Figure 15: Influence of HCUSP dissipation coefficients on hypersonic flow over 2-D wedge. 
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Figure 16: Influence of free-stream Mach number on the inviscid flow around NAGA 0012 
airfoil with the preconditioned HCUSP scheme. 
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Figure 17: Influence of preconditioning on the HCUSP scheme. 
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8 Concluding Remarks 

The CUSP scheme has been studied and analyzed. A detail comparison has been made 
between the CUSP , MATD, and scalar dissipation schemes. For transonic inviscid flows 
the CUSP scheme allows better resolution of shock waves, since they are captured with one 
interior point. However, the aerodynamic quantities such as lift and drag obtained with the 
original CUSP scheme are not as accurate on coarser meshes (i.e., 320 X 64 cells or less) 
as those calculated with the MATD scheme. Both the CUSP and MATD formulations can 
give high accuracy in the computation of viscous flows. In the case of high Re number flow 
over a flat plate, each of these schemes required only eight points in the boundary layer to 
have errors in computed skin-friction, displacement thickness, and momentum thickness that 
do not exceed 3%. Four times as many points is necessary to obtain comparable accuracy 
with the scalar scheme. For transonic viscous flows and coarser meshes the accuracy in 
aerodynamic coefficients is somewhat better with the MATD scheme than with the original 
CUSP scheme. This loss in accuracy with the CUSP scheme on coarser grids appears to be 
a consequence of the limiter producing a first-order scheme over significant portions of the 
flow held and higher levels of background dissipation. 

Modifications to the CUSP scheme for improving the coarse-grid accuracy have been 
presented. These changes restrict the activation regions of the first-order scheme to the 
neighborhoods of shock waves according to (3.2) and reduce background dissipation using 
the limiter of (3.8). They allow the CUSP scheme to give comparable accuracy to that 
obtained with the MATD scheme on coarse meshes. With these modifications to the CUSP 
scheme, convergence stall has been removed. Convergence has been further improved by 
introducing the aspect-ratio scaling factor of (2.24). 

In comparison to the scalar scheme the CUSP scheme requires roughly 25% more com- 
puter time while the MATD scheme needs about 15% more time. In general, convergence 
behavior with the CUSP and MATD schemes is similar. 

With our present choice of HCUSP dissipation coefficients it has been shown that the 
resolution of strong shock waves occurring in hypersonic flows is possible whereas the sim- 
plified coefficients that were published previously failed. At this point the HCUSP scheme 
appears to be a better choice than the present MATD scheme for hypersonic flow problems. 
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